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Abstract 

Biological  cells  are  composed  of  many  subsystems  and  organelles.  The  subsystem  called  the 
cytoskeleton  is  composed  of  long  rod-shaped  filaments.  They  give  the  cell  form  and  help  attach 
the  cell  to  the  substrate  and  neighbors.  One  type  of  the  filaments  is  made  up  of  a protein  called 
actin.  In  studying  cells  biologists  use  microscopes  that  can  be  focused  at  different  levels  and  can 
be  automated  to  take  multiple  images.  Various  stain  treatments  are  used  to  bring  out  various 
cell  characteristics.  In  this  study  a computational  method  is  studied  that  automatically  isolates 
the  actin  structures  in  cells.  The  algorithm  begins  by  segmenting  the  cell  with  a scheme  called  a 
watershed.  Due  to  noise  and  a large  number  of  local  minima  the  number  of  resulting  segmented 
regions  can  be  large  and  not  informative.  The  object  then  is  to  merge  regions,  based  on  nearness 
and  common  properties,  into  regions  called  features.  The  merging  approach  used  here  involves  a 
graph  representation  called  an  adjacency  graph  and  a search  algorithm  that  finds  the  connected 
components  in  the  graph.  These  connected  components  become  the  merged  features. 

Keywords:  actin  fibers;  biological  cells;  depth- first  search:  merging;  undirected  graph: 
watershed. 

1 Introduction 

Biological  cells  are  composed  of  many  subsystems,  called  organelles.  The  subsystem  of  most  interest 
in  this  work  is  called  the  cytoskeleton,  which  is  composed  of  long  rod-shaped  molecules,  or  filaments, 
that,  are  attached  to  each  other  and  to  other  organelles.  They  give  the  cell  form  and  help  attach 
the  cell  to  the  substrate  and  neighbors.  The  cytoskeleton  consists  of  three  types  of  filaments:  actin 
filaments,  intermediate  filaments,  and  microtubules.  For  a more  complete  discussion  of  the  cell 
cytoskeleton  and  its  structure  see  Bao  and  Suresli  [2]  and  Ethier  and  Simmons  [3],  pp.  23-29. 
Graphically,  actin  filaments  can  be  visually  identified  in  specially  stained  cell  images.  However, 
to  automatically  identify  cell  features,  such  as  actin  filaments,  is  not  straightforward.  Modern 
microscopic  systems  can  scan  more  cell  samples  than  individuals  can  visually  examine  to  compare,  for 

‘Applied  and  Computational  Mathematics  Division,  Information  Technology  Laboratory,  National  Institute  of 
Standards  and  Technology,  Gaithersburg  MD  20899-8910,  david.gilsinn@nist.gov. 

'Biochemical  Science  Division,  Material  Measurement  Laboratory,  National  Institute  of  Standards  and  Technology, 
Gaithersburg  MD  20899-8313.  kiran.bhadriraju@nist.gov. 

+ Biochemical  Science  Division.  Material  Measurement  Laboratory.  National  Institute  of  Standards  and  Technology, 
Gaithersburg  MD  20899-8313.  john.elliott@nist.gov 


1 


example,  cell  lines  so  it  is  important  that  computational  methods  be  developed  that  locate  features 
like  actin  filaments  so  that  changes  done  to  actin  shape  within  the  cell  edges  can  be  automatically 
identified. 

This  work  develops  an  initial  approach  to  locating  features  in  a cell  image.  It  combines  a method 
of  image  segmentation  into  what  will  be  called  regions,  with  an  algorithm  to  merge  regions  into 
more  defineable  cell  features.  The  algorithm  developed  is  not  the  only  approach  to  cell  feature 
identification,  but  appears  to  be  an  initial  workable  solution  to  automatic  cell  feature  identification. 
Further  work  is  clearly  needed  in  the  future. 

The  feature  identification  process  involves  a multi-stage  algorithm.  The  basic  idea  is  to  first  apply 
an  image  segmentation  algorithm,  in  the  current  study  called  watershedding,  that  provides  as  output 
a matrix  the  size  of  the  image,  called  a label  matrix.  The  matrix  elements  are  uniquely  numbered 
so  that  all  the  associated  image  pixels  in  a region  are  given  a unique  number  identifying  the  region. 
Thus,  for  example,  all  of  the  pixels  in  region  3 are  identified  with  3 in  the  associated  elements  of  the 
label  matrix.  The  segmentation  procedure  usually  generates  a large  number  of  regions  because  of 
all  of  the  local  minima  and  noise  in  the  image.  The  idea  behind  the  second  stage  of  the  algorithm  is 
to  merge  regions  with  similar  properties  in  order  to  isolate  features.  To  do  this,  a measure  is  defined 
that  is  used  to  link  regions  with  similar  properties.  In  mathematical  graph  theory  this  measure  is 
used  to  form  what  is  called  an  adjacency  graph  where  the  graph  nodes  are  the  numbered  regions  and 
the  link  weights  are  the  measures  of  adjacency.  For  a basic  introduction  to  graph  theory  see  Ore 
[8].  Once  an  adjacency  graph  is  developed  a tolerance  is  specified  that  tells  which  graph  links  are 
critical.  It  leads  to  a reduced  graph.  Not  all  nodes  are  connected  by  a path  in  the  graph  to  every 
other  node.  Thus,  given  a node,  a search  of  the  graph  is  made  for  all  the  nodes  linked  by  a path  to 
that  node.  These  node  groups  are  called  connected  components  and  they  ultimately  are  identified 
as  the  features.  Once  the  connected  components  of  the  graph  are  identified,  the  node  regions  can 
then  be  merged  by  combining  region  pixels.  Finally,  a new  label  marix  is  developed  that  is  used 
by  a boundary  tracing  algorthm  to  identify  the  feature  boundaries.  This  feature  label  matrix  could 
also  be  used  to  assist  in  post- processing  feature  analysis.  This  report,  however,  will  concentrate  on 
a description  of  the  implementation  of  the  merging  algorithm  in  MATLAB. 

Combining  watershed  segmentation  with  region  merging  is  not  a new  idea.  See.  for  example, 
Haris  et  al.  [6]  and  Pires  et  al.  [9].  What  is  missing  from  many  of  these  other  studies,  however, 
are  details  about  code  implementation.  The  current  study  is  an  attempt  to  fill  the  implementation 
gap  and  provide  a detailed  description  of  a code  implementation  using  MATLAB.  It  is  hoped  that 
this  will  provide  a basis  for  future  extensions  or  constructive  modifications  of  the  algorithm  and 
computer  code. 

It  is  assumed  that  the  reader  is  familar  with  basic  image  processing  terminology  and  is  proficient 
with  MATLAB.  A working  knowledge  of  the  MATLAB  Image  Processsing  Toolbox  and  some  graph 
theory  is  also  helpful. 

The  report  is  structured  as  follows.  Each  section  will  describe  in  detail  a portion  of  the  Waterslied- 
Merge  program  given  in  Appendix  A.  Section  2 presents  the  order  in  which  functions  are  called 
in  the  program.  Section  3 describes  the  image  contrast  adjustment  by  histogram  equalization, 
edge  sharpening  of  the  image  by  filtering,  and  finally  the  initial  segmentation  of  the  image  by 
watershedding  into  regions.  Section  4 describes  the  region  formation  in  terms  of  a data  type,  called 
a structure,  the  formation  of  the  region  boundary  lists,  creating  the  region  adjacency  structure, 
and  identifying  of  the  closest  adjacent  neighbors.  The  section  then  introduces  a graph  theoretic 
formation  of  the  adjacency  structure  that  simplifies  the  feature  and  feature  boundary  formation. 
Section  5 presents  images  of  the  results,  with  features  outlined  along  with  a table  of  merging  and 
timing  results.  Section  6 describes  areas  for  future  work.  Section  7 acknowledges  the  image  sources 
and  Section  8 is  the  formal  disclaimer  of  the  use  of  commercial  software. 
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2 Main  Watershed-Merge  Function 

This  function  organizes  the  overall  sequence  of  function  calls.  It  calls  each  of  the  functions  Pre.processor, 
Watershed.  Merge,  feature_boundary,  and  Post  .processor  in  turn.  The  term  function  in 
MATLAB  may  be  referred  to  in  other  programming  languages  as  subroutines.  Defining  the  main 
function  as  a function  type  allows  all  of  the  other  functions  to  be  included  in  the  same  hie  and  allows 
for  program  portability. 


3 Pre-processing  an  Image 

Segmentation  algorithms  can  tend  to  oversegment  an  image  into  a large  number  of  regions  that 
make  the  process  of  joining  similar  regions  to  form  a common  feature  difficult  and  time  consuming. 
Depending  on  the  number  of  oversegmented  regions  one  can  run  into  out-of- memory  errors  or  overrun 
allocated  arrays.  It  is  therefore  important  in  the  pre-processing  of  an  image  to  obtain  as  few  noise- 
based  regions  as  possible.  In  this  section  we  describe  a modified  watershed  approach  that  produces 
an  initial  set  of  segmented  regions  and  allows  a merging  algorithm  to  determine  potential  cell  features 
in  a reasonable  amount  of  compute  time. 

3.1  Histogram  Equalization 

An  image  can  be  considered  as  a large  matrix,  with  M rows  and  N columns,  where  the 

element  values,  or  pixels  for  picture  elements,  can  vary  from  0 to  255  for  8-bit  images  or  0 to  65535 
for  16-bit  images.  Let  nk.  k = 0,  ■ - •,  L — 1,  where  L is  255  for  8-bit  images  or  65535  for  16-bit 
images,  be  the  number  of  pixels  having  intensity  level  k.  If  the  range  of  intensity  values  is  divided 
into  intervals  or  ‘bins’  then  a plot  of  the  numbers  of  pixel  values  falling  within  the  bin  intervals  is 
called  a histogram.  It  shows  the  distribution  of  the  intensity  values  in  the  image.  In  particular,  a 
histogram  is  a discrete  function  h(rk)  = rik  where  rk  is  the  k- th  intensity  gray  scale  value  and  nk 
is  the  number  of  pixels  in  the  image  with  intensity  value  rk.  A normalized  histogram  is  given  by 
a function  p(rk)  = rik/MN  for  k — 0,  — 1.  p{rk)  can  be  thought  of  as  an  estimate  of  the 

probability  of  occurence  of  the  intensity  value  rk . 

The  distribution  of  pixel  intensities  in  an  image  can  vary  depending  on  the  image  acquisition 
system.  If  the  pixel  intensities  are  concentrated  in  the  lower  values  the  image  shows  up  as  ’dark', 
and  if  they  are  concentrated  in  the  upper  values  it  shows  up  as  ’light'.  The  object  of  histogram 
equalization  is  to  increase  the  intensity  range  as  much  as  possible  over  the  range  of  intensities. 
Figure  1 looks  totally  black  because  all  of  the  figure  intensities  are  concentrated  in  the  low  range. 
This  is  clear  from  the  histogram  of  the  intensities  shown  in  Fig.  2.  Notice  the  distribution  of 
the  intensities  in  the  dark  region  of  the  gray  scale  color  bar  at  the  bottom  of  the  figure.  This 
accounts  for  the  black  image.  The  image  is  intended  to  be  one  of  a biological  cell,  but,  due  to  lack 
of  contrast,  it  does  not  show  in  the  image.  The  initial  image  used  in  this  study  is  referenced  as 
celll_MLCPPandfactin_red_0075. 

Histogram  equalization  is  a process  that  modifies  the  histogram  of  an  image  in  order  to  achieve 
a better  overall  contrast.  It  is  shown  in  Gonzalez  and  Woods  [4],  pp.  122-128,  that  the  histogram 
equalization  of  the  histogram  intensities  can  be  given  by 

k L-  1 k 

sk  = T(rk)  = (L  - l)^p(rj)  = J2nJ-  k = 0.---.L  -1,  (1) 

j= o j= o 
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Figure  1:  Raw  Figure  without  Histogram  Equalization.  Intensities  are  Distributed  in  the  Low  Range. 
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Figure  2:  Histogram  of  the  Raw  Figure  without  Histogram  Equalization.  The  corresponding  inten- 
sities are  shown  as  gray  levels  at  the  bottom  of  the  figure. 
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Figure  3:  Raw  Figure  with  Histogram  Equalization.  Intensities  Redistributed  across  the  Range. 


where  n = no  + n\  + • ■ • 4-  ul- i is  the  total  number  of  pixels  and  Sk  is  the  redistributed  number  of 
pixels  with  intensities  equal  to  r in  the  bins.  The  default  number  of  bins  used  in  MATLAB  is  taken 
to  be  256.  Figure  3 shows  an  image  of  the  cell  after  histogram  equalization.  Note  the  redistribution 
of  the  pixel  intensities  towards  the  lighter  side  of  the  range  in  Fig.  4.  If  10  is  the  initial  image  then 
in  MATLAB  histogram  equalization  can  be  performed  by  the  call  II  = histeq(IO),  where  II  is  the 
resulting  image.  It  can  be  displayed  by  a call  to  imshow(Il)  if  desired. 

In  the  next  section  we  will  sharpen  the  fibers  a bit. 

3.2  Image  Filtering  and  Feature  Sharpening 

Although  histogram  equalization  leads  to  an  image  with  a wider  distribution  of  intensities  it  some- 
times introduces  some  blurring  of  edges.  There  is  definite  edge  blurring  of  the  cytoskeleton  fibers 
in  Fig.  3 as  a result  of  the  histogram  equalization.  Image  filters  can  be  applied  to  sharpen  the 
edges.  In  particular,  we  apply  a filter  called  an  unsharp  mask.  Unsharp  masking  is  a method  of 
edge  enhancement  in  an  image  that  involves  subtracting  a scaled  blurred  version  of  the  original 
image  from  the  original  image.  This  has  the  affect  of  enhancing  the  edges  in  the  image.  It  can  be 
implemented  by  subtracting  a 3 x 3 blurring  mask  from  a scaled  3x3  mask  with  a one  in  the  center 
and  zeros  elsewhere.  This  can  be  done  in  MATLAB  with  a call  to  the  function  fspecial  followed  by 
a call  to  imfilter.  It  is  illustrated  below,  where  h is  the  resulting  mask  and  0.2  is  the  scale  factor 
for  the  mask.  This  scale  factor  was  found  to  be  optimum  in  McAndrew  [7],  pp.  105-107.  The  call 
imfilter  applies  the  mask  to  the  image.  For  a brief  discussion  of  unsharp  masking  see  McAndrew 
[7],  The  image  boundary  is  also  padded  with  zeros  so  that  the  application  of  other  image  filters  will 
not  cause  errors  at  the  boundaries.  This  is  done  with  the  MATLAB  function  padarray. 

%unsharp  masking  is  edge  crisping  by  subtracting  a scaled  blurred  version 
°/0of  the  image  from  the  original 
h = f special (’ unsharp 0 . 2) ; 

12  = imfilter (II ,h) ; 
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Figure  4:  Histogram  of  the  Raw  Figure  with  Histogram  Equalization. 


°/0Display  the  sharpened  image 
figure,  imshow(I2) , impixelinf o ; 

°/0Pad  the  array  boundary  with  zeros 
12  = padarray(I2 , [1  1] ,0, ’both’ ) ; 

An  example  of  an  unsharp  filter  is  given  in  McAndrew  [7],  pp.  105-107  as 
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where  c is  a constant  chosen  to  provide  the  best  sharpening  result.  In  the  code  segment  above 
c = 0.2.  The  result  of  applying  the  unsharp  filter  to  the  image  in  Fig.  3 is  shown  in  Fig.  5.  The  cell 
fibers  have  been  sharpened  somewhat. 

3.3  Image  Segmentation  by  Watershedding 

In  order  to  understand  watershedding,  an  image  needs  to  be  thought  of  as  a surface  with  intensities 
as  surface  heights.  The  surface  will  have  peaks  and  valleys  as  a result  of  the  intensity  distribution. 
If  we  think  of  water  falling  on  the  surface  then  it  will  collect  in  valleys  around  the  local  minima  until 
it  overflows  into  neighboring  valleys.  The  overflows  occur  at  the  ridges  that  surround  the  valley 
around  a local  minimum.  The  watershed  transform  finds  the  local  valleys  and  marks  the  ridges  of 
the  valleys.  The  MATLAB  call  to  watershed  produces  a matrix  the  size  of  the  original  image  in 
which  each  matrix  element  is  associated  with  an  image  pixel.  The  matrix  produced  is  called  a label 
matrix  and  identifies  all  of  the  individual  valleys  around  minima  with  a unique  number  with  the 
elements  in  the  label  matrix  associated  with  the  pixels  in  each  valley  being  given  the  same  unique 
number.  The  ridges  between  valleys  are  labeled  with  zero.  For  a more  complete  discussion  of  the 
watershed  transform  see  Gonzalez  and  Woods  [4],  pp.  769-776. 
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We  illustrate  below  what  the  structure  of  a label  matrix  padded  with  zeros  might  look  like  for 
a simple  example  of  a 12  x 12  image  with  4 watershedded  regions.  Note  the  region  boundaries 
separated  by  zeros  and  the  outer  edges  padded  with  zeros.  We  will  use  this  example  throughout  the 
text  to  illustrate  some  steps  in  the  code. 
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The  one  difficulty  with  applying  the  watershed  transform  directly  to  an  image  is  that  it  can 
oversegment  an  image  into  too  many  regions  that  de-enhance  image  understanding  and  the  isola- 
tion of  significant  features.  This  is  illustrated  in  Fig.  6.  Although  the  segmentation  is  valid,  the 
oversegmentation  observed  in  Fig.  6 is  the  result  of  a large  number  of  local  minima  often  caused  by 
image  noise.  To  control  the  oversegmentation  we  first  apply  a smoothing  operation. 

We  will  briefly  discuss  a method  to  control  oversegmentation  described  in  Gonzalez  and  Woods 
[4],  pp.  776-778,  and  in  Gonzalez  et  al.  [5],  pp.  422-424.  Their  approach  to  minimizing  the 
number  of  valleys  around  minima  is  to  note  that  many  minima  are  shallow  relative  to  neighboring 
pixels  and  should  not  be  treated  as  legitimate  deep  valleys.  Such  minima  are  extraneous  to  the 
general  feature  identification.  The  method  described  in  Gonzalez  and  Woods  [4]  and  Gonzalez  et 
al.  [5]  is  based  on  a concept  of  markers.  Markers  are  connected  components  of  an  image.  There 
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Figure  6:  The  Watershed  Transform  applied  to  the  Original  Cell  Image. 


are  internal  markers  associated  with  objects  of  interest  and  external  markers  associated  with  the 
background.  An  effective  way  to  mimimize  the  effects  of  small  local  mimima  is  to  first-  apply  a 
smoothing  filter.  Once  the  image  has  been  smoothed  the  internal  markers  can  be  identified,  where 
an  internal  marker  has  these  characteristics  (1)  the  region  is  surrounded  by  higher  intensity  pixels, 
(2)  the  region  is  connected,  and  (3)  all  points  in  the  connected  region  have  the  same  intensity.  The 
internal  markers  can  be  identified  using  the  Image  Processing  Toolbox  (IPT)  in  MATLAB  by  a 
function  called  imextendedmin.  This  function  computes  the  set  of  minima  in  an  image  that  are 
lower  than  neighbors  by  a specified  threshold.  The  threshold,  in  a sense,  instructs  the  function  to 
ignore  all  minima  with  intensities  less  than  the  threshold  and  makes  the  internal  marker  identify  the 
locations  of  significant  local  minima.  The  process  is  called  an  extended  minima  transform. . For  a 
complete  discussion  of  extended  minima  see  Soille  [10],  pp.  203-204.  In  the  case  of  16-bit  cell  images 
with  intensity  range  of  0 to  65535,  a reasonable  choice  of  threshold  was  found  to  be  10000.  This 
threshold  was  found  by  applying  the  extended  minima  transform  to  multiple  images  and  noting  the 
reduction  in  the  initial  segmentation  of  the  act-in  fiber  regions. 

Once  the  internal  markers  have  been  identified,  the  watershed  algorithm  in  IPT  can  be  applied. 
This  is  done  in  Gonzalez  et  al.  [5],  pp.  422-425,  by  first  finding  the  pixels  in  the  image  generated 
by  the  imextendedmin  call  that  are  midway  between  each  of  the  internal  markers.  This  is  done 
by  applying  the  watershed  algorithm  to  an  image  that  is  the  distance  function  value  between  each 
pixel  in  the  internal  marker  image  and  the  nearest  nonzero  pixel.  In  the  case  below,  the  ’cityblock' 
metric  is  used,  which  is  the  sum  of  the  lengths  of  the  orthogonal  edges  of  the  right  triangle  formed 
between  two  pixel  coordinates.  Although  there  are  multiple  choices  of  metrics  for  the  MATLAB 
function  bwdist,  this  metric  seems  to  produce  good  results  for  the  cell  images. 

°/0Use  markers  to  eliminate  noise  and  shallow  minima  and  then  apply  watershed 
imin  = imextendedmin(I2 , 10000) ; 

L = watershed(bwdist(imin, ’cityblock’)) ; 

The  watershed  boundaries  here  are  the  external  markers.  The  resulting  overlays  to  the  original 
cell  image  can  be  displayed  by  the  function  imshow  in  IPT  by  giving  the  boundary  zeros  in  the 
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Figure  7:  The  Watershed  Transform  applied  to  the  Distance  Transform  of  the  Internal  Marker 
Image. 


label  matrix  a very  low  intensity  value.  The  segment  of  code  below  locates  the  label  matrix  elements 
that  are  zero.  It  then  copies  the  original  image  to  another  one  that  is  overlayed  with  low  intensity 
pixel  values  at  the  pixels  associated  with  zeros  in  the  label  matrix.  The  final  pre-processed  image  is 
thus  created  by  assigning  a very  low.  in  this  case  1.  intensity  value  to  the  watershed  boundary  pixels, 
identified  by  0 in  the  label  matrix  L.  For  the  original  cell,  this  is  given  in  Fig.  7.  This  process  has 
reduced  somewhat  the  segmented  regions  generated  by  applying  the  watershed  algorithm  directly 
to  the  original  image  and  makes  it  easier  to  merge  regions. 

'/(Display  the  region  boundaries  as  black  overlays  to  12 
em  = (L  ==  0)  ; 

13  = 12; 

13 (em)  = 1; 
figure,  imshow(I3); 
dravnow ; 


4 Region  Merging 

If  we  examine  Fig.  7 closely  we  can  note  that  many  of  the  regions  have  similar  properties.  In 
particular,  there  are  a number  of  regions  that  clearly  cover  areas  of  bundled  actin  filaments.  These 
filaments  surround  a middle  section  of  the  cell.  The  middle  section  of  the  cell  seems  to  have  similar 
properties  with  lower  intensities  and  not  many  filaments.  It  would  then  seem  that  we  would  like  to 
merge  regions  that  have  high  intensity  filament  structure,  and  also  merge  the  regions  covering  the 
middle  of  the  cell.  The  question  is  how  best  to  do  this. 

One  approach  that  has  been  used  in  the  literature  (Haris  et  al.  6 ) is  to  fink  neighboring  regions 
with  common  properties.  A standard  mathematical  structure  that  can  be  employed  as  a compu- 
tational tool  that  allows  efficient  representations  of  these  relations  is  called  a graph  (see  Ore  [8j. 
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Graphs  are  graphical  representations  that  join  points,  called  nodes , with  links  that  may  or  may  not 
have  some  form  of  weight  or  distance  associated  with  them.  Nodes  usually  have  numbers  associ- 
ated with  them.  Ordinarily,  links  with  large  weights  are  considered  the  most  potentially  significant 
candidates  for  merging.  An  example  might  be  a street  network  where  nodes  are  intersections  and 
streets  are  the  links  with  speed  limits  on  them.  Higher  speeds  mean  shorter  travel  times.  Another 
example  might  be  a communication  network  with  message  capacities  on  the  links.  Larger  capacities 
allow  for  higher  volume  message  passing.  In  this  study,  however,  we  will  take  links  with  low  weights 
to  be  significant.  That  is,  adjacent  regions  with  lower  weights  on  their  links  will  be  candidates  for 
merging. 

Graphs  are  convenient  tools  for  linking  information.  It  sometimes  happens  that  a complete 
description  of  a graph  may  not  be  necessary.  For  example,  two  nodes  with  similar  properties  may 
have  a link  between  them  with  a very  low  weight.  These  two  nodes  could  be  merged.  It  is  this  idea 
that  is  behind  the  process  of  merging  regions  in  a segmented  image.  Each  region  is  considered  to 
be  a node  and  links  are  established  between  neighboring  regions.  Regions  with  similar  properties 
are  considered  as  candidates  for  merging.  A graph  representation  of  regions  that  link  neighboring 
regions  is  called  an  adjacency  graph.  The  approach  to  region  merging  taken  in  this  work  combines 
an  adjacency  graph  representation  of  the  relationship  between  regions,  where  adjacency  is  defined  in 
terms  of  a distance  function  value  assigned  to  each  graph  link,  with  a boundary  detection  algorithm 
to  create  the  merged  region  boundaries. 

The  merging  process  begins  by  using  the  fact  that  the  watershed  algorithm  produces  a label 
matrix  that  assigns  to  each  located  region  a unique  number  and  identifies  the  boundary  of  the 
region  with  zeros  in  the  label  matrix  (see  (3)).  In  the  first  step  the  pixels  in  each  region  are  located 
along  with  the  boundaries  of  each  region  using  an  8-nrask,  a 3 x 3 matrix  centered  about  a pixel 
element  in  the  label  matrix.  Next,  neighboring  regions  are  located  by  scanning  around  the  boundary 
of  each  region  with  a 4-mask  in  order  not  to  asssociate  two  regions  as  neighbors  if  there  is  only  one 
common  boundary  point.  A 4-mask  is  represented  by  the  elements  above,  below  and  to  the  left  and 
right  sides  of  a pixel  element  in  the  label  matrix.  A distance  function  between  any  two  regions  is 
then  constructed  based  on  the  difference  between  the  average  intensities  of  the  two  regions.  Given 
that  region  neighbors  and  link  distances  are  known,  the  adjacency  graph  can  be  built  as  a matrix 
in  sparse  form  that  has,  for  each  row,  a region  and  a neighbor  along  with  the  link  distance.  Thus 
the  adjacency  graph  is  represented  in  node-node-link  weight  form.  The  adjacency  graph  is  in  an 
undirected  form  in  that  it  includes  links  from  A to  B and  B to  A,  as  well  as  links  A to  A.  Two  steps 
are  then  performed  to  compress  the  adjacency  graph.  The  first  is  to  eliminate  all  A to  A links  and 
the  second  is  to  eliminate  all  links  with  weights  greater  than  a predefined  tolerance.  The  adjacency 
matrix  is  then  searched  to  form  a linked  list  of  regions  to  be  merged.  Finally,  the  linked  list  is  used 
to  join  the  linked  region  pixels  into  what  is  called  a feature  and  a boundary  algorithm  is  used  to 
locate  the  external  boundary  of  the  features. 

4.1  Initial  Region  Structure  Formation 

In  order  to  understand  how  region  information  is  stored  one  needs  to  know  how  an  image  is  repre- 
sented as  a matrix.  The  origin  of  the  coordinate  system  used  to  locate  pixels  begins  in  the  upper 
left  corner  with  the  coordinate  (1,  1).  The  x-coordinate  points  to  the  right  and  is  indexed  by  matrix 
columns  and  the  y-coordinate  points  down  and  is  indexed  by  rows.  This  can  be  a bit  confusing  since 
in  numerical  matrix  analysis  one  thinks  in  terms  of  rows  first  then  columns,  such  as  (i,j).  But  in 
imaging  it  is  more  convenient  to  think  in  terms  of  columns  then  rows,  such  as  (j,i). 

Before  a region  data  structure  can  be  defined,  the  number  of  regions  found  by  watershed  can 
be  determined  by  the  call 


10 


n_reg  = max (max (L) ) ; 

where  L is  the  label  matrix  generated  from  the  watershed  call. 

The  data  representation  used  for  regions  in  this  study  is  called  a structure.  These  are  MATLAB 
arrays  with  attached  named  fields.  Structures  are  also  commonly  used  in  programming  languages, 
such  as  C.  Fields  can  contain  different  classes  of  data.  For  example,  one  held  might  contain  text, 
another  might  contain  a matrix.  MATLAB  allows  the  programmer  to  build  a structure  dynamically 
in  that  one  does  not  have  to  preallocate  the  structure  array  lengths.  To  begin  constructing  the  region 
structures  here,  a call  is  made  to  the  function  regionprops  in  the  IPT.  which  builds  a structure 
from  the  label  matrix  created  by  watershed.  In  the  call  below  we  want  to  store  as  fields  the  pixel 
list  of  a region  as  well  as  the  number  of  pixels  in  the  boundary  pixel  list.  Since  the  regionprops 
function  produces  a limited  number  of  fields,  we  will  construct  other  Helds  dynamically  as  we  need 
them. 

R = regionprops (L, ’PixelList ’ , 3 Perimeter’ ) ; 

If  i is  a region  index  then  the  initial  region  structure  looks  like 

R(i) . PixelList 
R(i) .Perimeter 

The  pixel  list  is  refered  to  in  (column,  row]  form  that  represents  an  (x.y)  coordinate  in  the  image 
framework. 

In  an  earlier  section  we  padded  the  original  image  with  zeros  around  the  boundary.  When 
watershed  is  applied  to  this  form  of  image  in  MATLAB  it  identifies  the  entire  image  boundary  as 
region  1.  We  ignore  region  1 by  overwriting  the  1‘s  in  the  label  matrix  with  zeros. 

col  = R(l) .PixelList (:, 1) ; 
row  = R(l) .PixelList (: ,2) ; 
for  rc  = 1 :length(col) 

L(row(rc) , col(rc))  = 0; 

end 

We  then  add  a pixel  intensity'  held  and  an  average  intensity'  held  to  the  region  structure.  In  order 
to  perform  arithmetic  on  pixel  elements  we  first  convert  the  image  to  a double  precision  matrix.  This 
code  segment  also  points  out  the  difference  between  image  pixel  coordinate  references  and  matrix 
element  references.  Note  that  the  c.r J pixel  reference  is  converted  to  an  r.c)  reference  in  the  matrix 
I2d  in  order  to  retrieve  the  matrix  element  for  the  pixel  intensity.  I2d  is  a matrix  the  size  of  the 
image  12  but  in  double  precision  form. 

/(Add  the  average  region  intensity  to  the  region  structure 

for  i = 2:n_reg 
sum  = 0.0; 

Pixel_num  = length (R(i)  .PixelListC:  , 1))  ; ‘/.get  the  number  of  pixels  in  R(i) 
for  j = l:Pixel_mm 

c = R(i) .PixelList (j ,1)  ; 
r = R(i) .PixelList (j ,21 ; 

R(i). Intensities (j)  = I2d(r,c); 

sum  = sum  + 10“  (-4)*I2d(r , zj ; /(scaling  16  bit  image  pixels 
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end 

R(i).avg_int  = (sum/Pixel_num) *10"4 ; 

end 

This  code  adds  to  the  region  structures  from  regions  2 to  n_reg  the  two  new  fields  called  Intensities 
and  avg-int.  The  region  structure  now  looks  like 

R(i) .PixelList 
R(i) . Perimeter 
R(i) . Intensities 
R(i) . avg_int 

In  computing  the  average  intensity  we  use  the  fact  that  intensity  values  are  of  the  order  of  104, 
in  a uintl6  image,  so  that,  if  they  are  scaled,  the  sums  can  be  formed  with  values  in  a much  lower 
range  and  then  scaled  back  after  the  scaled  average  has  been  computed.  This  is  needed  since  the 
number  of  pixels  in  a region  can  be  large.  The  sum  of  unsealed  numbers  could  be  very  large,  leading 
to  potential  numerical  inaccuracies  due  to  roundoff  or  overflow. 

4.2  Building  Region  Boundary  Lists 

To  build  the  region  boundary  lists,  we  first  have  to  estimate  an  upper  bound  for  all  of  the  boundary 
perimeters.  In  fact,  we  double  the  maximum  perimeter  length  in  order  to  make  sure  that  we  do  not 
overrnn  an  array.  This  length  is  used  in  a temporary  working  array  that  is  deleted  once  the  region 
boundary  lists  have  been  formed.  Note  that  we  make  sure  that  the  maximum  boundary  length  BL 
is  an  integer  by  using  the  MATLAB  fix  function. 

for  i = 2:n_reg 

Bdry_length(i)  = R(i) .Perimeter ; 

end 

"/Get  an  upper  bound  on  the  length  of  the  boundary  arrays 

BL  =2* (f ix(max (Bdry_length) ) +1) ; "/.possible  overestimate  but  is  safe 

We  can  now  add  a boundary  list  field  to  the  region  structure  by  scanning  the  region  pixels  with 
an  8-mask  and  locating  the  pixels  in  the  label  matrix  with  zero  values  neighboring  region  pixels. 
Points  in  a mask  are  referred  to  in  compass  notation:  W,  NW.  N,  NE,  E.  SE,  S,  SW.  If  a region 
pixel  value  is  designated  by  (r,c)  in  the  label  matrix  L.  the  NW  corner  of  an  8-mask  is  (r-l.c-1), 
which  uses  the  first  element  in  the  mask_offsets_8  array.  Similarly  the  N point  is  (r-l.c)  and  this  uses 
the  second  element  in  mask_offsets_8  array.  This  continues  around  the  mask.  We  scan  the  pixels 
in  each  region  from  2 to  n_reg.  As  we  scan,  we  may  pick  up  duplicate  boundary  array  values,  but 
these  are  eliminated  using  the  MATLAB  function  unique  that  eliminates  duplicate  values  in  a set. 
The  method  involved  with  the  mask  scanning  is  based  on  ideas  in  the  m-function  boundaries  given 
in  Gonzalez  et  al.  [5],  pp.  557-560.  The  script  below  sets  up  an  8-mask  pixel  offset  matrix.  The 
pixel  list  for  each  region  is  scanned.  A check  is  made  so  that  the  image  boundary  is  not  overrun. 
The  boundary  pixels  found  are  put  into  a temporary  array  that  is  used  in  the  function  unique  to 
eliminate  duplicate  pixels.  Finally  a boundary  list  is  formed  and  added  to  the  region  structure.  At 
the  end,  the  boundary  list  is  combined  with  the  region  pixel  list  to  form  a final  region  pixel  list. 


°/0Set  up  a temporary  working  array 
Temp  = zeros (BL, 2); 

'/.Next  array  to  hold  the  lengths  of  the  region  boundaries 
B = zeros (n_reg, 1) ; 

'/.Set  up  offsets  from  a pixel  in  (r,c)  form  to  locate  8-mask  entries  around 
°/0the  (r,c)  pixel. 

mask_of f sets_8  = [-1  -1;  -1  0;  -1  +1;  0 +1;  +1  +1;... 

+1  0;  +1  -1;  0 -1];  ‘/.for  boundary  detect  mask 

'/.Scan  each  region  pixel  to  locate  boundary  pixels  that  are  identified  as 
'/.zeros  in  the  label  matrix. 

for  i = 2:n_reg 

Pixel_num  = length(R(i)  .PixelList( : , 1))  ; '/.get  the  number  of  pixels  in  R(i) 
work_ count  = 0; 
for  j = l:Pixel_num 

c = R(i) . PixelList  (j , 1)  ; 
r = R(i) .PixelList (j ,2) ; 

'/.scan  with  8-point  mask  to  get  all  0’s  around  region 
'/.scan  NW  N NE  E SE  S SW  W 

for  k = 1:8  ’/.need  to  check  that  indices  are  not  out  of  range 

if  (r  + mask_off sets_8(k, 1)  <=  0)  I (c  + mask_of f sets_8(k, 2)  <=  0) I . . . 
(r  + mask_of f sets_8(k,  1)  > mL)  I (c  + mask_of f sets_8 (k , 2)  > nL) 
continue ; 

end 

if  L(r  + mask_off sets_8(k, 1) , c + mask_off sets_8(k,2))==  0 
‘/.This  will  likely  generate  duplicate  pixels  due  to 
'/.overlapping  scan  masks 
work_count  = work_count  + 1; 

Temp(work_count ,1)  = c + mask_off sets_8(k,2) ; 

Temp(work_count ,2)  = r + mask_off sets_8(k, 1) ; 

end 

end 

end 

C = unique  (Temp , ’ rows  ’)  ; '/.Eliminates  duplicate  boundary  points 

'/.Add  a boundary  list  to  the  structure  of  R 

B(i,l)  = length(C(  : , 1) ) ; '/.Boundary  length  for  R(i) 

R(i) .BoundaryList  = zeros (B (i) , 2) ; 

R(i)  .BoundaryList  (1  :B(i)  , 1)  = C(l:B(i),l);  '/.column 

R(i)  . BoundaryList  (l:B(i)  ,2)  = C(l:B(i),2);  '/.row 

clear  Temp  C 

end 
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7„Join  pixels  and  boundaries  to  essentially  eliminate  the  zeros  in 
0/„the  label  matrix.  The  only  zeros  left  are  the  image  boundaries. 
°/0Union  also  eliminates  duplicates. 


for  i = 2:n_reg 

R(i)  . PixelList  = union(R(i)  . PixelList  , R(i)  . BoundaryList , ’ rows  ’ ) ; 

end 


At  this  point  we  have  extended  the  region  structure  by  adding  a boundary  list. 


R(i) .PixelList 
R(i) . Perimeter 
R(i) . Intensities 
R(i) . avg_int 
R(i) .BoundaryList 


We  note  that  the  union  of  the  boundary  pixels  with  region  pixels  would  have  the  effect  of 
eliminating  the  boundary  zeros  as  shown  in  the  matrix  example  below. 


0 0 0 
0 3 3 
0 3 3 
0 3 3 
0 3 3 
0 3 4 

0 4 4 

0 4 4 

0 4 4 

0 4 4 

0 1 1 
0 0 0 


0 0 0 
3 3 2 

3 3 2 
3 3 2 

3 3 4 

4 4 4 

4 4 4 
4 4 4 
4 4 4 
4 4 4 

1 1 1 
0 0 0 


0 0 0 
2 2 2 
2 2 2 
2 2 2 
4 4 2 
4 4 4 
4 4 4 
4 4 4 
4 4 4 
1 1 1 
1 1 1 
0 0 0 


0 0 0 
2 2 0 
2 2 0 
2 2 0 
2 2 0 
2 2 0 
2 2 0 
4 2 0 
1 1 0 
1 1 0 
1 1 0 
0 0 0 


(4) 


4.3  Building  a Region  Adjacency  Structure 

We  now  have  a boundary  list  for  each  region,  although  we  ignore  region  1 since  the  watershed 
call  makes  region  1 the  image  boundary.  To  form  a neighbor  structure  for  a region  we  trace  the 
boundary  of  each  region  with  a 4-mask  and  check  the  region  numbers  above,  below,  and  to  the  sides 
of  the  mask.  The  scan  procedure  used  is  similar  to  the  one  used  to  find  region  boundaries,  but  in 
this  case  we  use  a 4-mask  to  identify  adjacent  regions.  This  mask  eliminates  neighboring  regions 
that  have  one  point  of  contact.  As  we  scan  around  a region  boundary  we  build  an  adjacency  array 
and  from  it  compute  the  number  of  neighbors  and  add  a neighbor  field  to  the  region  structure.  We 
note  that  the  array  B ( i ) is  computed  in  the  previous  script  portion. 


°/0At  this  point  we  have  a boundary  list  and  boundary  length  for  each  region 
%We  ignore  region  1 since  the  watershed  makes  that  the  image  boundary 
°/oNext  we  need  to  establish  an  Adjacency  Matrix  but  in  sparse  form. 

7oTo  do  this  we  trace  the  boundary  of  each  object  with  a 4-mask  and  check 
%the  numbers  above,  below  and  to  the  sides  of  the  mask. 
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test  = zeros (4,1); 


%Set  up  mask_of f set_4  for  N E S W Use  a 4-mask  in  order  not  to  link  regions 
%with  only  one  boundary  point  in  common 

mask_of f sets_4  = [0  -1;  1 0;  0 1;  -1  0];  %Mask  to  determine  neighbors 
for  i = 2:n_reg 

Adjacency (n_reg, 1)  = 0;  %Set  up  an  adjacency  array  for  a region 

for  j = 1 : B ( i)  °/0step  around  the  length  of  the  boundary  for  R(i) 

c = R(i) . BoundaryList ( j , 1) ; 
r = R(i) . BoundaryList (j ,2) ; 
for  k = 1:4 

rk  = r + mask_of f sets_4(k,  1)  ; 
ck  = c + mask_of f sets_4(k, 2) ; 

if  ((rk  ~=  0)  & (ck  ==  0))  I ((rk  ==  0)  & (ck  ~=0))... 

I ((rk  ==  0)  & (ck  ==0))|  rk  > mL  I ck  > nL 
continue;  '/.when  at  a boundary  element  or  out  of  range 

else 

test(k,l)  = L(rk,ck); 
for  m = 2:n_reg 

if  test(k,l)  ==  m;  %may  hit  m more  than  once. 

Adjacency (m, 1)  = 1; 

end 

end 

end 

end 

end 

°/0We  fill  in  the  R(i)  structure  elements  for  numbers  of  neighbors 
sum_neigh  = 0; 
for  isum  = l:n_reg 

sum_neigh  = sum_neigh  + Adjacency (isum, 1)  ; 

end 

number_neighbors  = sum_neigh; 

R(i) ,no_neighbors  = number_neighbors ; 

count  = 0; 

for  ii  = 2:n_reg 

if  (Adjacency (ii , 1)  ~=  0) 
count  = count  + 1 ; 

R(i) .neighbor (count)  = ii; 

end 

end 

clear  Adjacency 

end 


15 


We  have  now  extended  the  structure  to  include  the  number  of  neighbors  and  the  neighbor  list. 

R(i) . PixelList 
R(i) . Perimeter 
R(i) . Intensities 
R(i) . avg_int 
R(i) .BoundaryList 
R(i) . no_neighbors 
R(i) .neighbor 

4.4  Identifying  Similar  Neighbors 

Although  we  have  a list  of  all  neighbors  of  a region,  not  all  neighbors  have  similar  characteristics.  In 
order  to  merge  regions  into  features  with  more  unique  characteristics,  we  need  to  define  a measure 
of  closeness.  For  this  we  define  a distance  function.  d(i,j)  between  regions  i and  j,  that  satisfy  the 
standard  properties  of  a metric. 

d(i, j)  >=  0, 

d(i,j)  = 0 if  and  only  if  i = j, 

d(i  , j ) = d( j , i) , 

d(i,j)  <=  d(i,h)  + d(h,j). 

The  distance  measure,  D(p,l),  we  use  in  this  study  involves  differences  between  the  average 
intensities  of  two  regions’  pixels.  Since  the  image  intensities  are  of  the  order  104  we  scale  quantities 
so  that  they  fall  in  the  interval  [0,1].  The  absolute  values  of  the  average  intensity  differences  are 
used  since  the  squares  of  the  difference  could  create  values  of  the  order  108.  We  first  extend  the 
region  structure  by  adding  a distance  value  between  each  region  i and  all  neighboring  regions, 
R(j).D  Jieighbor(i).  Then  we  define  a tolerance,  Tol.  that  is  a constraint  on  the  measure  of  the 
closest  neighbors  in  terms  of  region  properties  and  extend  the  region  structure  by  adding  the  closest 
neighbors  to  the  structure.  A similar  measure  has  been  used  by  Haris  et.  al.  [6].  Pires  et  al.  [9] 
add  a term  involving  a boundary  gradient.  Other  distance  functions  could  be  defined  based  on  such 
things  as  region  texture  measures,  but  the  one  used  here  seems  to  provide  reasonable  results  and 
has  shown  success  in  the  literature. 

Once  the  distance  function  has  been  defined,  the  fields  for  the  number  of  adjacent  regions  within 
a distance  tolerance  are  computed.  The  scaled  tolerance  used  involves  the  median  scaled  by  1/1.1, 
This  scale  was  determined  experimentally  and  appears  to  produce  the  most  meaningful  feature 
formation.  An  adjacency  matrix  is  formed  with  number  of  rows  equal  to  the  number  of  regions  and 
column  width  the  size  of  the  list  of  maximum  number  of  adjacent  regions  within  tolerance  to  some 
region.  This  matrix  is  sparse  and  is  used  in  a search  algorithm  to  find  all  connected  components  in  the 
adjacency  graph.  The  neighbors  within  tolerance  to  region  j are  stored  in  R(j).neighJn_Tol(p)  and 
the  number  of  the  neighbors  within  tolerance  are  stored  in  R(j).no_neigh _in_Tol.  This  information 
is  used  to  create  the  adjacency  matrix,  adj. 

70The  next  step  is  to  convert  the  region  structure  to  a generalized 
"/Adjacency  matrix  in  sparse  form  using  a distance  measure  between  regions 
"/This  is  where  intensities  enter  the  picture.  This  section 
°/0of  code  produces  a matrix  of  paths  from  an  origin  nose  to  a final 
"/destination . Region  merging  will  be  performed  by  averaging  values  for  the 
"/region  features  for  the  regions  in  each  path. 
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°/0We  will  use  a process  of  origin-destination  link  matrices  to  create 
°/„a  final  graph 

p = 0;  “/(need  to  accumulate  over  all  regions 
for  j = 2:n_reg 

n_neigh  = R( j ) . no_neighbors ; 
for  i = l:n_neigh 
P = p+i; 

J(p,l)  = j;  %get  column 

I(p,l)  = R(j)  .neighbor(i)  ; “Aget  row 

“/(compute  distance  function  as  10“  (-4)  times 

“/(average  intensity  differences  since  image  intensities  are 

“/(order  10~4 

D(p,l)  = (10~(-4))*abs((R(j) .avg_int  - R(I (p) ) . avg_int) ) ; 

R(j) ,D_neighbor(i)  = D(p,l); 

end 

end 

°/„At  this  point  we  need  to  estimate  the  merging  tolerance.  Experimentation 
“/(suggests  the  following  is  a working  tolerance 
Tol  = median(D) /I . 1 ; 

“/(Add  to  the  structure  the  number  of  neighbors  within  tolerance 
for  j = 2:n_reg 

n_neigh  = R( j ) . no_neighbors ; 

P = 0; 

for  i = l:n_neigh 

RD  = R( j ) . D_neighbor (i) ; 
if  (RD  <=  Tol)&(RD  > eps) 
p = p+i; 

R( j ) . neigh_in_Tol (p)  = R(j) .neighbor(i) ; 

end 

end 

R( j ) . no_neigh_in_Tol  = p; 

end 

“/0Make  an  adjacency  list  matrix  in  a modified  linked  form.  Each  row 
°/0will  begin  with  a region  and  the  subsequent  columns  will  list  all  of  the 
“/(regions  within  the  distance  tolerance  of  that  region.  This  matrix  will  be 
%used  in  a search  algorithm  during  the  merge  phase.  It  is  global. 

no_adj  = zeros(l ,n_reg) ; 
no_adj(l)  = 0; 
ad j (1,1)  = 0; 
for  i = 2:n_reg 

no_adj(i)  = R(i) . no_neigh_in_Tol ; 
for  j = l:no_adj(i) 

adj(i,j)  = R(i) .neigh_in_Tol(j) ; 

end 
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end 


4.5  Building  an  Undirected  Graph  Representation  of  the  Region  Adja- 
cency 

We  can  now  represent  the  adjacency  between  regions  using  an  undirected  graph  representation  in 
terms  of  nodes  and  links.  In  the  present  case  nodes  will  represent  regions  and  links  will  be  identified 
in  terms  of  two  nodes.  Each  link  has  a distance  measure  between  nodes,  or  regions.  The  undirected 
graph  is  represented  in  terms  of  a sparse  matrix  where  each  row  has  the  form  of  node-node-link 
weight,  with  the  distance  measure  for  the  link  weight.  We  then  cull  the  matrix  of  all  links  of  a node 
to  itself  and  links  with  distances  greater  than  the  specified  tolerance.  Finally  a list  of  all  of  the 
critical  nodes  is  formed. 

%We  create  an  undirected  graph  of  neighbors. 
length_I  = length(I); 

%find  all  node  pairs  of  the  same  value,  such  as  1 1 or  2 2,  etc. 

kill  = zeros ( 1 , length_I) ; 

Index  = find(D  < eps) ; 
kill(Index)  = 1; 
clear  Index 

y„Form  a reduced  undirected  node-node-link  weight  matrix  by  killing  node 
°/„pairs  of  equal  value  index 

p = 0; 

for  i = l:length_I 
if  (kill (i)  ~=  1) 

P = p+i; 

GDCp.l)  = j(i); 

GD (p , 2)  = I (i)  ; 

GD (p , 3)  = D(i)  ; 

end 

end 

%Next  eliminate  links  with  link  weights  greater  than  Tol. 

[ng,mg]  = size(GD); 
q = 0; 

for  j = 1 : ng 

if (GD( j ,3)  < Tol) 
q = q+1; 

GDT(q, 1)  = GD(j ,1) ; 

GDT(q , 2)  = GD(j  ,2) ; 

GDT(q,3)  = GD ( j , 3) ; 

end 

end 
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clear  GD 


"/All  of  the  critical  nodes  in  the  adjacency  graph  will  appear  in  GDT(:,1). 

’/.All  other  nodes  in  the  set  difference  between  2:ndeg  and  GDT(:,1)  will  be 
"/.taken  as  separate  features.  The  critical  nodes,  called  V,  will  be  used  in  a 
"/subsequent  search  algorithm  during  the  region  merge  phase. 

V = unique  (GDT( 1)  ) ; "/.get  the  unique  nodes. 

4.6  Feature  Formation 

The  vector  V may  not  contain  the  indices  of  all  regions.  We  identify  those  regions  not  appearing 
in  V as  single  region  features.  Thus  features  contain  one  or  more  regions.  In  forming  features  we 
will  begin  by  creating  the  single  region  features  and  then  create  the  multiple  region  features.  In 
the  process  of  creating  features  we  apply  a graph  search  algorithm,  DFS,  for  depth-first  search. 
In  doing  so,  we  need  to  create  an  array  that  tells  us  whether  we  have  visited  a prior  region  in  the 
search  process  so  that  we  don’t  end  up  in  a never  terminating  graph  loop.  This  is  the  array  visited. 
The  array  TotaLvisited  lists  all  of  the  nodes  visited.  We  update  this  array  as  we  search  through 
the  graph  for  connected  components.  The  connected  components  are  identified  as  the  features. 
Connected  components  are  formed  by  starting  at  a node  that  has  not  been  visited  and  going  to  each 
of  its  linked  nodes  in  the  adj  matrix  and  then  following  down  the  graph  from  each  of  those  nodes. 
It  is  clear  that  this  becomes  a recursive  process  that  ends  when  there  are  no  more  linked  nodes  that 
have  not  been  visited.  The  search  algorithm  is  implemented  in  the  depth-first  search  function  DFS. 
Once  the  connected  components  for  a node  have  been  identified,  a feature  structure  can  be  formed. 

The  depth-first  search  algorithm  is  based  on  the  algorithm  SEARCH  in  Aho  et  al.  [1].  p.  176- 
179.  It  takes  as  global  variables  visited,  nomdj,  adj.  V and  begins  with  a node  called  start  and 
sets  visited(start)=  1.  The  start  node  is  always  taken  as  the  lowest  index  numbered  region,  say  r, 
such  that  TotaLvisited  = 0.  Then  it  loops  through  each  node  adjacent  to  start.  Let  one  be  called 
w.  If  visited(w)  = 1 it  goes  on  to  another  adjacent  node  to  start.  If  not  then  it  calls  DFS(w) 
recursively.  In  effect  it  works  down  all  path  links  for  the  adjacent  node  w.  When  all  the  path  linked 
nodes  to  w have  been  visited,  it  proceeds  to  the  next  adjacent  node  to  start.  Again  it  continues 
until  all  possible  path  linked  nodes  to  start  have  been  visited.  One  of  the  significant  outputs  of  this 
portion  of  the  code  is  a new  label  matrix  that  labels  the  features  in  the  same  fashion  that  the  region 
label  matrix  did  for  the  regions.  That  is,  all  of  the  pixels  in  each  feature  are  uniquely  numbered. 

"/At  this  point  we  have  an  undirected  graph  with  adjacency  distances  all 
"/within  tolerance.  We  will  begin  the  feature  merging  by  marking  those 
"/regions  that  have  been  visited  during  the  merging  process 

Total_visited  = zeros (1 ,n_reg) ; 

"/Find  all  regions  that  do  not  appear  in  GDT(:,1).  First  set 
°/Total_visited(l)  = 1 to  eliminate  it  from  being  a region 

Total_visited(l)  = 1; 

°/We  can  find  all  of  the  single  node  region  features  by  a set  difference 
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Total_nodes  = 2:n_reg; 

Single_reg  = setdif f (Total_nodes , V) ; 

‘/.Form  features  for  each  of  the  Single-nodes 

ISn  = length(Single_reg) ; 

‘/.Make  the  first  ISn  features  the  same  as  the  regions  identified  in  the 
°/.Single_nodes  list.  We  form  a feature  structure  that  contains  the  feature 
‘/.pixels . 

for  i = l:lSn 

Sn  = Single_reg(i) ; 

lplist  = length(R(Sn)  . PixelList  ( : , 1) ) ; °/0get  no  pixels  in  reg  Sn 
‘/.build  the  feature  structure 

F(i) . PixelList (1 : lplist , 1)  = R(Sn) . PixelList (1 : lplist , 1) ; 

F(i) . PixelList (1 : lplist , 2)  = R(Sn) .PixelList(l  :lplist ,2) ; 
Total_visited(Sn)  = 1; 

end 

°/0Initialize  counter  for  the  next  multi-node  features 
p = ISn; 

°/0Initialize  a new  label  Matrix 
NewL  = zeros (mL ,nL) ; 

‘/.Before  creating  features  we  need  to  find  the  regions  that  are  to  be  merged 
‘/.for  each  feature 

I = f ind(Total_visited  ==  0);  ‘/.Find  all  nodes  that  haven’t  been  visited 

while  'isempty(I)  ‘/.Continue  as  long  as  there  are  nodes  that  haven’t  been 
‘/.visited 

visited  = zeros ( 1 ,n_reg) ; 

start  = 1(1);  ‘/.Use  the  first  available  0 visited  node 

‘/.The  next  call  does  a depth-first  search  of  the  undirected  graph.  The 

‘/.output  is  the  array  visited  of  all  nodes  ultimately  linked  to  the  node 

‘/.start . This  generates  the  connected  components  to  start  in  the  graph. 

DFS(start);  ‘/.get  the  connected  components  linked  to  node  start 

clear  I ‘/.I  will  be  updated  below 

p = p+1;  ‘/.start  a new  feature. 

‘/.Initialize  the  new  feature 

Reg_visit  = find(visited  ==  1);  ‘/.Find  the  regions  visited  in  the  search 
IRv  = length(Reg_visit) ; 

F(p) . PixelList  = R(Reg_visit ( 1) ). PixelList ; 
for  i = 1 : IRv  ‘/.Join  all  the  pixels  of  visited  regions 
F(p) . PixelList  = union(F(p) . PixelList ,.. . 
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R(Reg_visit (i) ) . PixelList , 'rows ’ ) ; 

end 

Total_visited  = Total_visited  + visited;  °/0add  recent  visited  to 

°/0total  visited 

I = f ind(Total_visited  ==  0); 

end 

no_features  = p 

%set  the  pixel  values  in  F(q) . PixelList  in  the  array  NewL. 

for  q = 1 : no_f eatures 

"/.form  a new  label  matrix 
new_col  = F(q) . PixelList (:,  1) ; 
new_row  = F(q) . PixelList (:, 2) ; 
for  nrc  = 1 : length(new_col) 

if  (new_row(nrc)  ==  0)  I (new_col (nrc)  ==  0) 
continue ; 

end 

NewL(new_row(nrc) ,new_col (nrc) ) = q; 

end 

end 

We  now  have  the  merged  pixel  list  for  each  feature  and  a new  label  matrix  for  all  of  the  features. 
The  feature  structure  at  this  point  consists  of  a pixel  list. 

F(i) . PixelList 

We  next  need  to  determine  the  boundaries  of  the  merged  regions.  For  this  we  apply  a boundary 
tracking  algorithm  by  Moore  (see  Gonzalez  and  Woods  [4],  p.  796-797). 

4.7  Feature  Boundary  Identification 

In  order  to  locate  the  features  on  the  original  image,  being  able  to  identify  the  feature  boundaries  is 
critical.  This  code  portion  describes  the  implementation  of  the  Moore  boundary  tracking  algorithm 
used.  The  basic  algorithm  is  rather  simple.  It  starts  in  the  upper  left  corner  of  a feature  and  sets  that 
pixel  as  the  first  point  of  the  boundary.  It  places  an  8-mask  over  the  point  and  begins  a clockwise 
test  of  each  element  of  the  mask  beginning  with  western  point  of  the  mask  until  the  test  hits  a pixel 
of  the  feature.  This  point  is  put  into  the  boundary  list  and  an  8-mask  is  set  over  this  point  and  the 
pixel  just  before  the  hit  would  be  outside  of  the  feature,  so  a clockwise  test,  beginning  at  the  last 
outside  point,  of  the  elements  of  the  8-mask  are  made  until  a pixel  hit  of  the  feature  is  made.  The 
process  continues.  The  stopping  criteria  is  a return  to  the  initial  pixel  and  a retest  around  the  initial 
pixel.  If  the  next  hit  is  the  same  as  the  hit  after  beginning  the  algorithm,  then  the  algorithm  stops, 
otherwise  it  continues.  There  are  cases  in  which  the  initial  pixel  is  hit  and  the  next  hit  is  not  the 
second  pixel  hit.  The  algorithm  adds  this  new  hit  to  the  boundary  list  and  continues.  For  a more 
complete  discussion  of  the  algorithm  see  Gonzalez  and  Woods  [4],  p.  796-798. 

The  final  output  of  this  portion  of  the  code  is  a new  label  matrix  with  feature  boundaries  identified 
by  zeros  in  the  same  manner  that  the  watershed  segmentation  algorithm  identified  boundaries  of 
the  segmented  regions.  This  label  matrix  is  used  to  plot  the  final  image  in  the  post-process  portion 
of  the  code. 
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"/  Start  the  Moore  Boundary  Algorithm 
[Nr,Nc]  = size(NewL); 
empty  = 0; 

ignore  = zeros (no_features , 1) ; 
for  p = 1 : no_f eatures 

"/Finding  the  feature  starting  boundary  point. 

°/0In  Matlab  arrays  are  stored  by  columns.  Thus  to  find  when  NewL  ==  p 
°/0the  result  returned  is  the  set  of  all  elements  in  linear  order 
7„starting  with  the  upper  left  point  of  the  feature.  This  point 
"/.represents  the  initial  boundary  point  for  the  feature. 

L_index  = min(f ind(NewL  ==  p));  '/Get  the  minimum  of  the  linear  indices 
if  isempty (L_index)  "/If  the  index  list  is  empty  get  a new  feature 
empty  = empty  + 1 ; 
ignore (empty)  = p; 
continue ; 

end 

"/Find  the  row  and  column  in  NewL  for  the  minimum  index 
r = mod(L_index,  Nr); 
if  r ==  0 
r = Nr; 

c = L_index/Nr; 

else 

c = f ix(L_index/Nr)  + 1; 

end 

°/0We  initialize  the  feature  boundary  structure 
F(p) .BoundaryList(l , 1)  = c; 

F(p) . BoundaryList (1 , 2)  = r; 
b_cnt  = 1 ; 

°/0By  selection  of  (r,c)  the  western  point  (r,c-l)  is  not  in  the  feature 
70or  on  the  boundary  of  the  pth  feature.  We  test  around  (r,c)  beginning 
°/0at  (r,c-l)  to  find  the  next  point  that  intersects  the  boundary.  If  no 
"/intersection  is  found  the  feature  is  a single  point.  The  function 
°/0Neigh8  generates  the  row  and  column  offset  given  a mask  index.  It  is 
°/0a  look-up  table, 
for  index  = 2:9 

[of f set_r , of f set_c]  = Neigh8(index) ; 
if  (NewL(r+offset_r,c+offset_c)  ~=  p) 
continue ; 

elseif  (NewL(r+off set_r ,c+offset_c)  ==  p) 

b_cnt  = 2;  '/Store  the  second  boundary  hit  for  stopping 
F(p) . BoundaryList (2 , 1)  = c+offset_c; 

F(p) . BoundaryList (2 , 2)  = r+offset_r; 

“/Get  the  point  before  the  hit 

[back_of f _r ,back_of f_c]  = Neigh8(index-1) ; 

Last_non_boundary  = [r+back_of f _r  , c+  back_off_c]  ; 
r = r+offset_r; 
c = c+offset_c; 
break;  "/exit  the  search  loop 
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end 


end 

"/If  the  search  loop  exits  with  index  = 9 then  it  has  returned  to  the 
°/„initial  non-boundary  point  and  this  feature  is  a one  point  feature, 
if  (index  ==  9) 

continue;  "/go  to  the  next  feature 

end 

"/To  continue  locating  the  boundary  we  reset  (r+of f set_r , c+of f set_c)  to 
°/„(r,c)  and  start  the  scan  around  (r,c)  from  the  Last_non_boundary 
°/„point.  We  set  a check  flag  to  be  used  as  a feature  boundary 
/(termination  check. 
check_cnt  =-999; 

while  (F(p) . BoundaryList (b_cnt , 1)  ~=  F(p) . BoundaryList (2 , 1) ) I . . . 

(F(p) . BoundaryList (b_cnt ,2)  ~=  F(p) . BoundaryList (2 , 2) ) | .. . 

(b_cnt  ~=  check_cnt  + 1)  "/.Termination  criteria 
"/search  around  the  boundary  point  center  (r,c)  beginning  with  the 
"/last  non-boundary  point.  The  function  cal  inv_Neigh8  takes  a row 
land  column  offset  and  produces  the  mask  start  index.  It  is  a 
°/look-up  table . 

offset_r  = Last_non_boundary (1)  - r; 
offset_c  = Last_non_boundary (2)  - c; 
start_index  = inv_Neigh8(offset_r ,offset_c) ; 
for  index  = start_index+l : start_index+8 

[of f set_r , of f set_c]  = Neigh8(index)  ; "/get  offsets  from  (r,c) 
if  (NewL(r+off set_r  ,c+of fset_c)  ~=  p) 
continue ; 

elseif  (NewL(r+of f set_r , c+of f set_c)  ==  p) 
b_cnt  = b_cnt  + 1; 

F(p) . BoundaryList (b_cnt , 1)  = c+offset_c; 

F(p) . BoundaryList (b_cnt ,2)  = r+offset_r; 

if  (F(p) .BoundaryList (b_cnt , 1)  ==  F(p) . BoundaryList ( 1 , 1) ) &.. 
(F(p) . BoundaryList (b_cnt , 2)  ==  F(p) . BoundaryList ( 1 , 2) ) 
check_cnt  = b_cnt; 

end 

"/Get  the  point  before  the  hit 

[back_of f _r ,back_of f _c]  = Neigh8(index-1) ; 

Last_non_boundary  = [r+back_of f _r , c+  back_off_c]  ; 

"/set  up  next  boundary  search  center 
r = F(p) . BoundaryList (b_cnt , 2) ; 
c = F(p) . BoundaryList (b_cnt , 1) ; 
break;  "/exit  the  local  boundary  search  loop 

end 

end 

end 

"/Delete  duplicate  points 

F(p) . BoundaryList  = unique (F(p) . BoundaryList rows ’ ) ; 

end 

"/To  display  the  merged  regions  first  form  a label  matrix  with  zeros  at 
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Figure  8:  Merged  Regions  into  Feat  ures  for  Cell  1. 


%the  feature  boundaries 
NewL_bdr  = zeros (Nr , Nc) ; 

NewL_bdr  = NewL; 
for  p = 1 : no_f eatures 

ign  = find (ignore  ==  p) ; 
if  "isempty (ign) 
continue ; 

end 

lFpB  = length(F(p) . BoundaryList ( : , 1) ) ; 
f or  q = 1 : lFpB 

r = F(p) .BoundaryList (q, 2) ; 
c = F(p) .BoundaryList (q,l) ; 

NewL_bdr(r ,c)  = 0; 

end 

end 

bdrs  = (NewL_bdr  ==  0) ; 

I2new  = 12; 

I2new(bdrs)  = 1 ; 
f igure , imshow(I2new) ; 

5 Results 

Figure  8 shows  the  reduced  number  of  regions  compressed  into  regions  that  outline  sections  of  the 
cell  with  common  feat  ures.  One  of  the  outlined  regions  is  the  cent  ral  port  ion  of  t he  cell  t hat  probably 
contins  the  nucleus.  There  are  other  regions  that  outline  the  dense  act  in  liber  bundles.  I here  are 
also  regions  that  identify  where  the  fibers  spread  out,  and  also  regions  that  outline  the  boundan 
portion  of  the  cell.  The  merging  algorithm  reduced  410  watershed  segmented  regions  to  1 It). 
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Cell.  Feature  Region  Nos. 

and  Timings  Is ) 

Cell 

Nr 

n2 

Tx 

t2 

t3 

t4 

T, 

1 

416 

110 

0.29 

0.73 

31.92 

1.10 

0.04 

2 

1024 

268 

0.11 

0.72 

42.11 

1.79 

0.04 

4 

2009 

435 

0.12 

0.73 

71.89 

2.71 

0.04 

5 

319 

75 

0.12 

0.72 

48.37 

0.59 

0.04 

6 

1155 

287 

0.12 

0.78 

43.79 

1.37 

0.04 

7 

1373 

343 

0.11 

0.74 

27.35 

1.60 

0.04 

8 

2193 

515 

0.11 

0.73 

47.63 

2.51 

0.04 

Table  1:  Cell  Regions.  Features,  and  Program  Timings 


Figure  b:  Raw  Image  for  Cell  2. 


We  display  below  the  image  results  of  actin  feature  identification  for  several  cells.  We  only  in- 
clude the  raw  cell  images  and  the  final  feature  identification  images.  The  following  are  the  reference 
file  names  for  the  cell  images:  cel!2_N ILCP Pandfacti  n_red_0078 . cell4WlLCPPand:acCn_red_0084. 
ceU5_MLCPPan  lfactin_red_0098.  cell6_MLC  PPandfaetin-recLi  1 01  .._JMLCPPandfactin_red_0104. 
cell8_MLCPPandfactin_red_0107. 

Table  1 has  8 columns  representing  cell  number,  watershed  regions,  feature  regions,  and  particular 
program  timings  in  seconds.  Column  1 Is  the  cell  number  associated  with  the  hie  number.  Column 
2.  is  the  number  of  watershed  segmented  regions.  Column  3.  IC-  Is  the  number  of  final  tea*  ire 
regions  after  merging.  Column  4.  Ti.  is  pre-process  timing.  Column  5.  T2.  is  Watershed  Crning. 
Column  6.  T*.  is  the  merge  process  timing.  Column  7.  T4.  is  the  feature  boundary  formation  by 
the  Moore  algorithm.  Column  8.  T5.  is  the  post-process  timing.  It  is  clear  that  future  work  must 
be  aimed  at  timing  reduction  in  the  merging  portion  of  the  code. 
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Figure  10:  Feature  Image  for  Cell  2. 


Figure  11:  Raw  Image  for  Cell  4. 
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Figure  12:  Feature  Image  for  Cell  4. 


Figure  13:  Raw  Image  for  Cell  5. 
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Figure  14:  Feature  Image  for  Cell  5. 


Figure  15:  Raw  Image  for  Cell  6. 
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Figure  16:  Feature  Image  for  Cell  6. 


Figure  17:  Raw  Image  for  Cell  7. 
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Figure  18:  Feature  Image  for  Cell  7. 


Figure  19:  Raw  Image  for  Cell  8. 
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Figure  20:  Feature  Image  for  Cell  8. 

6 Future  Work 

The  current  code  is  intended  as  a benchmark  in  that  its  results  can  be  compared  against  those  of 
other  feature  identification  algorithms.  There  are  several  directions  that  future  work  can  take.  First, 
there  are  some  code  inefficiencies  that  need  to  be  adjusted.  Other  segmentation  algorithms  could  be 
implemented  but  their  outputs  would  have  to  include  a region  label  matrix.  Other  distance  measures 
for  region  adjacency  could  be  experimented  with  as  well  as  other  distance  tolerances.  A multiple 
pass  through  the  algorithm  might  further  feature  compression.  The  current  code  is  structured  as  a 
one-pass  algorithm.  One  other  possible  efficiency  might  be  to  use  a modified  version  of  the  Moore 
boundary  tracking  algorithm  in  the  process  of  constructing  the  initial  region  boundaries. 
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A Watershed_Merge  MATLAB  Program 

function  Main_Watershed_Merge 

0/„MAIN_WATERSHED_MERGE  - This  function  simply  calls  the  functions  that 
'/.create  image  regions,  merge  them  into  features,  and  displays  the  final 
'/.image 

y. 

'/.Author : 

% David  E.  Gilsinn 

% Mathematical  and  Computational  Sciences  Division 
% National  Institute  of  Standards  and  Technology 
7.  100  Bureau  Drive,  Stop  8910 

“/.  Gaithersburg,  MD  20899-8910 

7. 

12  = Pre_processor ; '/.Get  the  padded  filtered  cell  image 

[L,I2]  = Watershed(I2)  ; °/0Get  the  label  matrix  of  the  region  boundaries 

[NewL,F,  no_features]  = Merge(L,  12);  '/.Merge  the  regions  into  features 
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NewL_bdr  = f eature_boundary (F , NewL,  no_features)  ; 7Locate  the  feature  boundaries 
Postprocessor (12,  NewL_bdr)  ; 7Plot  the  feature  boundaries  as  image  overlays 

%************************************************************************* 

7Main  functions 

************************************************************************ 
function  12  = Pre_processor 

7PRE_PR0CESS0R  - Loads  a cell  image,  sharpens  the  edges  and  pads  the 
7 boundary  with  zeros 

7 

70utput : 

7 12  - filtered  image  with  padded  boundary 

7 Plotted  image 

7 

7Author : 

7 David  E.  Gilsinn 

7 Mathematical  and  Computational  Sciences  Division 
7 National  Institute  of  Standards  and  Technology 

7 100  Bureau  Drive,  Stop  8910 

7 Gaithersburg,  MD  20899-8910 

7 

7Load  a 16  bit  image 

710  = imread( ’ celll_MLCPPandf actin_red_0075 . tif ’ ) ; 716  bit  image 
710  = imread( ’ cell2_MLCPPandf actin_red_0078 . tif ’ ) ; 716  bit  image 
710  = imread( ’ cell3_MLCPPandf actin_red_0081 . tif ’ ) ; 716  bit  image 
710  = imread( ’ cell4_MLCPPandf actin_red_0084 . tif ’ ) ; 716  bit  image 
710  = imread( ’ cell5_MLCPPandf actin_red_0098 . tif ’ ) ; 716  bit  image 
710  = imread( ’ cell6_MLCPPandf actin_red_0101 . tif ’ ) ; 716  bit  image 
710  = imread( ’ cell7_MLCPPandf actin_red_0104 . tif ’ ) ; 716  bit  image 

10  = imread( ’ cell8_MLCPPandf actin_red_0107 . tif ’ ) ; 716  bit  image 
7Do  a histogram  equalization 

11  = histeq(IO) ; 

7unsharp  masking  is  edge  crisping  by  subtracting  a scaled  blurred  version 
7of  the  image  from  the  original  image 
h = fspecial( ’unsharp’ ,0 . 2) ; 

12  = imf ilter (II ,h) ; 

7Display  the  sharpened  image 
figure,  imshow(I2) , impixelinf o ; 

7Pad  the  array  boundary  with  zeros 
12  = padarray(I2,  [1  1] , 0 , ’ both ’ ) ; 

function  [L , I 2]  = Watershed(I2) 

7WATERSHED  - This  function  makes  use  of  watershed  segmentation  of  an  image 
7into  regions.  The  function  initially  applies  an  image  extended 
7minimization . See  the  Soille  reference. 

7 
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°Z Input : 

7c  12  - Image  padded  with  zeros  around  the  image  boundary 
7c 

“/Output : 

7c  12  - The  loaded  image,  assumed  to  be  uintl6 

7c  L - label  matrix  the  size  of  12  with  watershed  regions  uniquely 
7c  identified  with  numbers.  That  is,  each  matrix  element  associated 

7c  with  a region  pixel  is  given  the  same  number.  Region  boundaries  are 

7c  identified  with  zeros 

7c 

“/.Author : 

7c  David  E.  Gilsinn 

7c  Mathematical  and  Computational  Sciences  Division 

7c  National  Institute  of  Standards  and  Technology 
7c  100  Bureau  Drive,  Stop  8910 

7c  Gaithersburg,  MD  20899-8910 

1 

“/Use  markers  to  eliminate  noise  and  shallow  minima  and  then  apply  watershed 
imin  = imextendedmin(I2 , 10000) ; 

L = watershed(bwdist (imin , ’ cityblock ’ ) ) ; 

“/Display  the  region  boundaries  as  black  overlays  to  12 
em  = (L  ==  0) ; 

13  = 12; 

13 (em)  = 1; 
figure,  imshow(I3); 
drawnow ; 

%************************************************************************** 
function  [NewL,F,  no_features]  = Merge(L,  12) 

“/MERGE  - This  function  takes  a label  matrix  of  a zero  padded  image  12, 
“/develops  a region  structure  array,  finds  the  region  boundaries  and  region 
“/neighbors  within  a specified  tolerance  distance,  forms  an  adjacency  array, 
“/and  links  neighbors  along  paths  from  a given  region  into  feature  regions. 

I 

“/Input : 

7c  12  - Zero  padded  16-bit  image 

7c  L - Label  matrix  of  regions  identified  with  unique  numbers 

7c 

“/Output : 

7c  NewL  - Label  matrix  for  the  features 

7c  F - Feature  structure  with  associated  fields 

7c  no_features  - number  of  features 

7c 

“/Author : 

7c  David  E.  Gilsinn 

7c  Mathematical  and  Computational  Sciences  Division 

7c  National  Institute  of  Standards  and  Technology 
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'/.  100  Bureau  Drive,  Stop  8910 

"/.  Gaithersburg,  MD  20899-8910 

1 

global  visited  no_adj  adj  V "/.Used  in  depth-first  search  algorithm 
[mL,nL]  = size(L); 

I2d  = double(I2);  "/.Convert  to  double  precision  for  computational  reasons 
n_reg  = max (max (L))  "/.Get  the  number  of  regions 

°/0Build  the  first  fields  of  the  region  structure 

"/.Begin  building  the  region  structure  array 

"/.Get  a list  of  the  pixels  and  approximate  boundary  length  of  the  objects 
"/. PixelList  is  given  in  (c,r)  form.  Need  to  zero  out  region  1. 

°/0Watershed  makes  the  zeros  around  the  boundary  of  the  image  region  1 . 

R = regionprops (L PixelList Perimeter ’) ; "/.R  structure 

'/.Zero  out  region  1,  although  this  might  not  be  necessary 

col  = R(l) . PixelList (:, 1) ; 
row  = R(l) . PixelList (:, 2) ; 
for  rc  = 1 : length(col) 

L(row(rc) ,col(rc))  = 0; 

end 

"/.We  only  work  with  region  2 to  n_reg. 

"/.Next  find  an  upper  bound  estimate  of  the  number  of  boundary  pixels 
for  i = 2:n_reg 

Bdry_length(i)  = R(i) . Perimeter ; 

end 

"/.Get  an  upper  bound  on  the  length  of  the  boundary  arrays 

BL  =2*(f ix(max(Bdry_length))+l)  ; "/.possible  overestimate  but  is  safe 

"/.Add  the  intensities  and  average  region  intensity  to  the  region  structure 
"/.Scale  the  16  bit  intensities  so  that  sume  can  be  formed  in  a reasonable 
"/.numerical  range.  16  bit  intensities  are  of  the  order  1CT4 

for  i = 2:n_reg 
sum  = 0.0; 

Pixel_num  = length(R(i)  .PixelListC  : ,1))  ; "/.get  the  number  of  pixels  in  R(i) 
for  j = l:Pixel_num 

c = R(i) . PixelList (j , 1) ; 
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r = R(i) . PixelList (j , 2) ; 

R( i) . Intensities (j ) = I2d(r,c); 

sum  = sum  + 10“ (-4) *I2d(r , c) ; ^scaling  16  bit  image  pixels 

end 

R(i).avg_int  = (sum/Pixel_num) *10“4 ; 

end 

70  ************************************************************************* 

°/0Find  the  region  boundary  and  add  the  field  to  the  region  structure 
7************************************************************************* 

7„Set  up  a temporary  working  array 

Temp  = zeros (BL, 2); 

7oNext  array  to  hold  the  lengths  of  the  region  boundaries 
B = zeros (n_reg , 1) ; 

7oSet  up  offsets  from  a pixel  in  (r,c)  form  to  locate  8-mask  entries  around 
%the  (r,c)  pixel. 

mask_of f sets_8  = [-1  -1;  -1  0;  -1  +1;  0 +1;  +1  +1;... 

+ 1 0;  +1  -1;  0 -1];  7.for  boundary  detect  mask 

7»Scan  each  region  pixel  to  locate  boundary  pixels  that  are  identified  as 
7oZeros  in  the  label  matrix. 

for  i = 2:n_reg 

Pixel_num  = length(R(i)  .PixelList( : , 1))  ; 7»get  the  number  of  pixels  in  R(i) 

work_count  = 0; 

for  j = l:Pixel_num 

c = R(i) . PixelList (j , 1) ; 
r = R(i) . PixelList (j , 2) ; 

7oScan  with  8-point  mask  to  get  all  0’s  around  region 
%scan  NW  N NE  E SE  S SW  W 

for  k = 1:8  70need  to  check  that  indices  are  not  out  of  range 

if  (r  + mask_offsets_8(k,l)  <=  0)  | (c  + mask_of f sets_8(k , 2)  <=  0)1. 
(r  + mask_of f sets_8(k, 1)  > mL)  I (c  + mask_off sets_8(k,2)  > nL) 
continue ; 

end 

if  L(r  + mask_off sets_8(k, 1) ,c  + mask_of f sets_8(k,2) )==  0 
7»This  will  likely  generate  duplicate  pixels  due  to 
7»overlapping  scan  masks 
work_count  = work_count  + 1; 

Temp(work_count , 1)  = c + mask_off sets_8(k,2) ; 

Temp(work_count ,2)  = r + mask_of f sets_8(k , 1) ; 


end 


end 


end 

C = unique  (Temp, ’rows’ ) ; "/.Eliminates  duplicate  boundary  points 

°/„Add  a boundary  list  to  the  structure  of  R 

B(i,l)  = length  (C  ( 1)  ) ; "/.Boundary  length  for  R(i) 

R(i) . BoundaryList  = zeros(B(i)  ,2)  ; 

R(i)  . BoundaryList ( 1 : B(i)  , 1)  = C(l:B(i),l);  "/.column 
R(i)  .BoundaryList  (1  :B(i)  ,2)  = C(l:B(i),2);  "/.row 

clear  Temp  C "/.clear  working  arrays  so  they  can  be  reset  for  each  region 

end 

"/.Join  pixels  and  boundaries  to  essentially  eliminate  the  zeros  in 
"/.the  label  matrix.  The  only  zeros  left  are  the  image  boundaries. 

"/.Union  also  eliminates  duplicates. 

for  i = 2:n_reg 

R(i) . PixelList  = union (R(i) . PixelList ,R(i) . BoundaryList , ’ rows ’ ) ; 

end 

yj**************************************************3*:*****!*:********1)'****!*:** 

"/.Build  the  sparse  form  of  the  adjacency  matrix 

"/.At  this  point  we  have  a boundary  list  and  boundary  length  for  each  region 
"/.We  ignore  region  1 since  the  watershed  makes  that  the  image  boundary 
"/.Next  we  need  to  establish  an  Adjacency  Matrix  but  in  sparse  form. 

"/.To  do  this  we  trace  the  boundary  of  each  object  with  a 4-mask  and  check 
"/.the  numbers  above,  below  and  to  the  sides  of  the  mask. 

test  = zeros(4,l); 

"/.Set  up  mask_of f set_4  for  N E S W Use  a 4-mask  in  order  not  to  link  regions 
"/.with  only  one  boundary  point  in  common 

mask_of f sets_4  = [0  -1;  1 0;  0 1;  -1  0];  “/.Mask  to  determine  neighbors 
for  i = 2:n_reg 

Adjacency (n_reg , 1)  = 0;  "/.Set  up  an  adjacency  array  for  a region 

for  j = l:B(i)  "/.step  around  the  length  of  the  boundary  for  R(i) 

c = R(i) . BoundaryList (j , 1) ; 
r = R(i) . BoundaryList (j ,2) ; 
for  k = 1:4 

rk  = r + mask_of f sets_4(k, 1) ; 
ck  = c + mask_off sets_4(k,2) ; 

if  ((rk  ~=  0)  & (ck  ==  0))  I ((rk  ==  0)  & (ck  ~=0))... 
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I ((rk  ==  0)  & (ck  ==0))|  rk  > mL  I ck  > nL 
continue;  ’/.when  at  a boundary  element  or  out  of  range 

else 

test(k,l)  = L(rk,ck); 
for  m = 2:n_reg 

if  test(k,l)  ==  m;  70may  hit  m more  than  once. 
Adjacency(m, 1)  = 1; 

end 

end 

end 

end 

end 

“/We  fill  in  the  R(i)  structure  elements  for  numbers  of  neighbors 
sum_neigh  = 0; 
for  isum  = l:n_reg 

sum_neigh  = sum_neigh  + Adj acency (isum , 1) ; 

end 

number_neighbors  = sum_neigh; 

R(i) . no_neighbors  = number_neighbors ; 

count  = 0; 

for  ii  = 2:n_reg 

if  (Adjacency (ii , 1)  ~=  0) 
count  = count  + 1 ; 

R(i) .neighbor (count)  = ii; 

end 

end 

clear  Adjacency  “/Array  is  reset  for  a neu  region 

end 

“/The  next  step  is  to  convert  the  region  structure  to  a generalized 
“/.Adjacency  matrix  in  sparse  form  using  a distance  measure  between  regions 
“/This  is  where  intensities  enter  the  picture.  This  section 
“/of  code  produces  a matrix  of  paths  from  an  origin  nose  to  a final 
“/destination.  Region  merging  will  be  performed  by  averaging  values  for  the 
“/region  features  for  the  regions  in  each  path. 

“/„We  will  use  a process  of  origin-destination  link  matrices  to  create 
“/a  final  graph 

p = 0;  “/need  to  accumulate  over  all  regions 
for  j = 2:n_reg 

n_neigh  = R( j ) . no_neighbors ; 
for  i = l:n_neigh 
P = P+i; 

J(p,l)  = j;  “/get  column 

I(p,l)  = R(j)  .neighbor(i)  ; “/get  row 

“/compute  distance  function  as  10~(-4)  times 

“/average  intensity  differences  since  image  intensities  are 
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7„order  10~4 

D(p,l)  = ( 10“ (-4) ) *abs ( (R ( j ) . avg_int  - R(I (p) ) . avg_int) ) ; 

R(j) .D_neighbor(i)  = D(p,l); 

end 

end 

°/„At  this  point  we  need  to  estimate  the  merging  tolerance.  Experimentation 
"/.suggests  the  following  is  a working  tolerance 
Tol  = median(D) /I . 1 ; 

"/.Add  to  the  structure  the  number  of  neighbors  within  tolerance 
for  j = 2:n_reg 

n_neigh  = R( j ) . no_neighbors ; 

p = 0; 

for  i = l:n_neigh 

RD  = R( j ) . D_neighbor (i) ; 
if  (RD  <=  Tol)&(RD  > eps) 

P = p+i; 

R(j) .neigh_in_Tol(p)  = R( j ). neighbor (i) ; 

end 

end 

R( j ) . no_neigh_in_Tol  = p; 

end 

"/.Make  an  adjacency  list  matrix  in  a modified  linked  form.  Each  row 
"/.will  begin  with  a region  and  the  subsequent  columns  will  list  all  of  the 
"/.regions  within  the  distance  tolerance  of  that  region.  This  matrix  will  be 
"/.used  in  a search  algorithm  during  the  merge  phase.  It  is  global. 

no_adj  = zeros (1 ,n_reg) ; 
no_adj (1)  = 0 ; 
ad j (1,1)  = 0; 
for  i = 2:n_reg 

no_adj(i)  = R(i) . no_neigh_in_Tol ; 
f or  j = 1 :no_adj (i) 

adj(i,j)  = R(i) . neigh_in_Tol (j ) ; 

end 

end 

%************************************************************************** 
"/.We  create  an  undirected  graph  of  neighbors  using  the  adjacency  matrix 

"/.We  create  an  undirected  graph  of  neighbors. 
length_I  = length(I) ; 

"/.find  all  node  pairs  of  the  same  value,  such  as  1 1 or  2 2 , etc. 
kill  = zeros (1 ,length_I) ; 
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Index  = find(D  < eps)  ; 
kill(Index)  = 1; 
clear  Index 

"/(Form  a reduced  undirected  node-node-link  weight  matrix  by  killing  node 
°/„pairs  of  equal  value  index 

P = 0; 

for  i = l:length_I 
if  (kill (i)  ~=  1) 
p = p+i; 

GD(p.i)  = J(i); 

GD (p , 2)  = I (i)  ; 

GD (p , 3)  = D(i) ; 

end 

end 

"/(Next  eliminate  links  with  link  weights  greater  than  Tol . 

[ng,mg]  = size(GD) ; 

q = 0; 

for  j = 1 : ng 

if (GDC j ,3)  < Tol) 
q = q+1; 

GDT (q , 1)  = GD(j  ,1)  ; 

GDT (q, 2)  = GD(j ,2) ; 

GDT (q, 3)  = GD( j ,3)  ; 

end 

end 

clear  GD 


°/oAll  of  the  critical  nodes  in  the  adjacency  graph  will  appear  in  GDT(:,1). 
"/(All  other  nodes  in  the  set  difference  between  2:ndeg  and  GDT(:,1)  will  be 
"/(taken  as  separate  features.  The  critical  nodes,  called  V,  will  be  used  in  a 
/(subsequent  search  algorithm  during  the  region  merge  phase.  It  is  global. 

V = unique  (GDT( 1) ) ; 70get  the  unique  nodes. 

°/„Build  the  merged  features  by  doing  a depth-first  search  of  the  unirected 
"/(graph 

%************************************************************************** 

°/„At  this  point  we  have  an  undirected  graph  with  adjacency  distances  all 
"/(within  tolerance.  We  will  begin  the  feature  merging  by  marking  those 
"/(regions  that  have  been  visited  during  the  merging  process 

Total_visited  = zeros(l ,n_reg) ; 
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%Find  all  regions  that  do  not  appear  in  GDT(:,1).  First  set 
%Total_visited(l)  = 1 to  eliminate  it  from  being  a region 

Total_visited(l)  = 1; 

%We  can  find  all  of  the  single  node  region  features  by  a set  difference 
7«between  all  nodes  and  V. 

Total_nodes  = 2:n_reg; 

Single_reg  = setdif f (Total_nodes , V) ; 

°/0Form  features  for  each  of  the  Single-nodes 

ISn  = length(Single_reg) ; 

%Make  the  first  ISn  features  the  same  as  the  regions  identified  in  the 
7oSingle_nodes  list.  We  form  a feature  structure  that  contains  the  feature 
%pixels . 

for  i = 1 : ISn 

Sn  = Single_reg(i) ; 

lplist  = length(R(Sn)  . PixelList  ( : , 1)  ) ; 7«get  no  pixels  in  reg  Sn 
7»build  the  feature  structure 

F(i) .PixelList (1 : lplist , 1)  = R(Sn) . PixelList ( 1 : lplist , 1) ; 

F(i) . PixelList ( 1 : lplist , 2)  = R(Sn) . PixelList ( 1 : lplist , 2) ; 
Total_visited(Sn)  = 1; 

end 

7»Initialize  counter  for  the  next  multi-node  features 
p = ISn; 

^Initialize  a new  label  Matrix 
NewL  = zeros(mL,nL) ; 

7oBefore  creating  features  we  need  to  find  the  regions  that  are  to  be  merged 
70for  each  feature 

I = f ind(Total_visited  ==  0);  7»Find  all  nodes  that  haven’t  been  visited 

while  ~isempty(I)  %Continue  as  long  as  there  are  nodes  that  haven’t  been 
70visited 

visited  = zeros(l ,n_reg) ; 

start  = 1(1);  %Use  the  first  available  0 visited  node 
7»The  next  call  does  a depth-first  search  of  the  undirected  graph.  The 
70output  is  the  array  visited  of  all  nodes  ultimately  linked  to  the  node 
7oStart . This  generates  the  connected  components  to  start  in  the  graph. 
7oThe  search  is  done  using  the  adjacency  matrix  adj  . 
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DFS (start);  °/0get  the  connected  components  linked  to  node  start 
clear  I °/J  will  be  updated  below 
p = p+1;  “/.start  a new  feature. 

^Initialize  the  new  feature. 

'/.The  visited  array  guards  against  cycles  in  the  graph  search. 

Reg_visit  = find(visited  ==  1);  °/0Find  the  regions  visited  in  the  search 
IRv  = length(Reg_visit) ; 

F(p) . PixelList  = R(Reg_visit ( 1) ) . PixelList ; 
for  i = 1 : IRv  '/.Join  all  the  pixels  of  visited  regions 
F(p) . PixelList  = union(F(p) . PixelList ,.. . 

R(Reg_visit (i) ) .PixelList , ’ rows ’ ) ; 

end 

Total_visited  = Total_visited  + visited;  '/.add  recent  visited  to 

'/.total  visited 

I = f ind(Total_visited  ==  0); 

end 

°/0Set  the  number  of  connected  components  as  the  number  of  features 
no_features  = p 

’/.set  the  pixel  values  in  F(q)  . PixelList  in  the  array  NewL. 

for  q = 1 : no_f eatures 

’/.form  a new  label  matrix 
new_col  = F(q) . PixelList (:, 1) ; 
new_row  = F (q) . PixelList (:, 2) ; 
for  nrc  = 1 : length(new_col) 

if  (new_row(nrc)  ==  0)  I (new_col (nrc)  ==  0) 
continue ; 

end 

NewL (new_row(nrc) ,new_col (nrc) ) = q; 

end 

end 


'/.Support  functions 

%************************************************************************* 
'/.Depth-first  graph  search  algorithm 

function  DFS(v) 

°/„DFS  - Depth-first-search  algorithm 

'/.See  Aho-Hopcroft-Ullman,  "The  Design  and  Analysis  of  Computer  Algorithms", 
°/.pp  176-177 

global  visited  no_adj  adj  V 
visited(v)  = 1; 
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na  = no_adj (v) ; 
for  i = l:na 

w = adj (v, i)  ; 
if  visited(w)  ==  0 
DFS(w) ; 

end 

end 

708-mask  offset  look-up  table 

function  [offset_r,  offset_c]  = Neigh8 (index) 

%NEIGH8  - Given  an  element  index  of  an  8-mask  this  function  produces  the 
%pixel  offset  from  the  central  pixel.  The  indices  are:  W - 1,  NW  - 2, 

7.N  - 3,  NE  - 4,  E - 5,  SE  - 6 , S - 7,  SW  - 8 
icase  = mod(index , 8) ; 
if  icase  ==  0 
icase  = 8; 

end 

switch  icase 
case  1 

offset_r  = 0; 
offset_c  = -1; 
case  2 

offset_r  = -1; 
offset_c  = -1; 
case  3 

offset_r  = -1; 
offset_c  = 0; 
case  4 

offset_r  = -1; 
offset_c  = 1; 
case  5 

offset_r  = 0; 
offset_c  = 1; 
case  6 

of f set_r  = 1 ; 
offset_c  = 1; 
case  7 

offset_r  = 1; 
offset_c  = 0; 
case  8 

offset_r  = 1; 
of f set_c  = -1 ; 

end 

%8-mask  index  look-up  table 

function  index  = inv_Neigh8(offset_r ,off set_c) 
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°/„INV_NEIGH8  - Given  a pixel  offset  an  index  is  produced 

“/.This  function  begins  with  neighborhood  offsets  and  produces  the  index 

“/.number  for  an  8 Neighbor  point.  The  indices  are:  W - 1,  NW-2,  N - 3, 


”/.NE  - 4,  E - 5, 
if  (offset_r  == 
index  = 1 ; 
elseif  (offset_r 
index  = 2; 
elseif  (offset_r 
index  = 3 ; 
elseif  (offset_r 
index  = 4 ; 
elseif  (offset_r 
index  = 5 ; 
elseif  (offset_r 
index  = 6 ; 
elseif  (offset_r 
index  = 7 ; 
elseif  (offset_r 
index  = 8; 


SE-6,  S - 7,  SW-8 
0)  & (offset_c  ==  -1) 

==  -1)  & (offset_c  ==  -1) 

==  -1)  & (offset_c  ==  0) 

==  -1)  & (offset_c  ==  1) 

==  0)  & (offset_c  ==  1) 

==  1)  & (offset_c  ==  1) 

==  1)  & (offset_c  ==  0) 

==  1)  & (offset_c  ==  -1) 


end 


function  NewL_bdr  = f eature_boundary(F,  NewL,  no_features) 

°/.FEATURE_BOUNDARY  - Form  the  feature  boundaries  and  add  a boundary  field 
“/.array.  The  algoritm  is  based  on  the  Moore  boundary  algorithm.  See 
“/.Gonzalez,  R.  C.,  Woods,  R.  E.,  Digital  Image  Processing,  Pearson-prentice 
“/.Hall,  Upper  Saddle  River,  NJ,  2008,  796-798. 

% 

“/.Input : 

“/.  F - Feature  structure 

“/.  NewL  - Feature  label  matrix 

“/.  no_features  - number  of  features 

% 

“/.Output : 

“/.  NewL-bdr  - Feature  label  matrix  with  boundaries  identified  by  zeros 

“/. 

“/.Author : 

“/.  David  E.  Gilsinn 

“/.  Mathematical  and  Computational  Sciences  Division 

“/.  National  Institute  of  Standards  and  Technology 
“/.  100  Bureau  Drive,  Stop  8910 

“/.  Gaithersburg,  MD  20899-8910 

“/. 

%************************************************************************** 

“/.  Start  the  Moore  Boundary  Algorithm 
[Nr,Nc]  = size(NewL); 
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empty  = 0 ; 

ignore  = zeros (no_features , 1) ; 
for  p = 1 : no_f eatures 

"/(Finding  the  feature  starting  boundary  point. 

"/(In  Matlab  arrays  are  stored  by  columns.  Thus  to  find  when  NewL  ==  p 
°/0the  result  returned  is  the  set  of  all  elements  in  linear  order 
°/„starting  with  the  upper  left  point  of  the  feature.  This  point 
"/(represents  the  initial  boundary  point  for  the  feature. 

L_index  = min(f ind(NewL  ==  p));°/0Get  the  minimum  of  the  linear  indices 
if  isempty (L_index)  °/0If  the  index  list  is  empty  get  a new  feature 
empty  = empty  + 1 ; 
ignore (empty)  = p; 
continue ; 

end 

"/(Find  the  row  and  column  in  NewL  for  the  minimum  index 
r = mod(L_index,  Nr); 
if  r ==  0 
r = Nr; 

c = L_index/Nr; 

else 

c = f ix(L_index/Nr)  + 1; 

end 

"/.We  initialize  the  feature  boundary  structure 
F(p) . BoundaryList (1 , 1)  = c; 

F(p) . BoundaryList (1 ,2)  = r; 
b_cnt  = 1; 

"/.By  selection  of  (r,c)  the  western  point  (r,c-l)  is  not  in  the  feature 
‘/,or  on  the  boundary  of  the  pth  feature.  We  test  around  (r,c)  beginning 
%at  (r,c-l)  to  find  the  next  point  that  intersects  the  boundary.  If  no 
".(intersection  is  found  the  feature  is  a single  point.  The  function 
°/0Neigh8  generates  the  row  and  column  offset  given  a mask  index.  It  is  a 
/(look-up  table . 
for  index  = 2:9 

[of f set_r , of f set_c]  = Neigh8 (index) ; 
if  (NewL(r+of f set_r , c+off set_c)  ~=  p) 
continue ; 

elseif  (NewL(r+offset_r ,c+offset_c)  ==  p) 
b_cnt  = 2; 

F(p) .BoundaryList (2,1)  = c+offset_c; 

F(p) .BoundaryList (2, 2)  = r+offset_r; 

°/0Get  the  point  before  the  hit 

[back_of f _r ,back_of f _c]  = Neigh8 (index-1) ; 

Last_non_boundary  = [r+back_of f _r , c+  back_off_c] ; 
r = r+offset_r; 
c = c+offset_c; 
break:;  "/(exit  the  search  loop 

end 

end 
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°/„If  the  search  loop  exits  with  index  = 9 then  it  has  returned  to  the 
/(initial  non-boundary  point  and  this  feature  is  a one  point  feature, 
if  (index  ==  9) 

continue;  °/0go  to  the  next  feature 

end 

/(To  continue  locating  the  boundary  we  reset  (r+of f set_r , c+of f set_c)  to 
%(r,c)  and  start  the  scan  around  (r,c)  from  the  Last_non_boundary 
/(point.  We  set  a check  flag  to  be  used  as  a feature  boundary 
/(termination  check. 
check_cnt  =-999; 

while  (F(p) . BoundaryList (b_cnt , 1)  ~=  F(p) . BoundaryList (2 , 1) ) I . . . 

(F(p) . BoundaryList (b_cnt , 2)  ~=  F(p) . BoundaryList (2 , 2) ) I .. . 

(b_cnt  “=  check_cnt  + 1) 

'/.search  around  the  boundary  point  center  (r,c)  beginning  with  the 

/(last  non-boundary  point 

offset_r  = Last_non_boundary (1)  - r; 

offset_c  = Last_non_boundary (2)  - c; 

start_index  = inv_Neigh8 (of f set_r ,of f set_c) ; 

for  index  = start_index+l : start_index+8 

[of  f set_r  , of  f set_c]  = Neigh8(index)  ; °/0get  offsets  from  (r,c) 
if  (NewL(r+of f set_r , c+of f set_c)  ~=  p) 
continue ; 

elseif  (NewL(r+of f set_r  , c+of f set_c)  ==  p) 
b_cnt  = b_cnt  + 1; 

F(p) . BoundaryList (b_cnt , 1)  = c+offset_c; 

F(p) . BoundaryList (b_cnt , 2)  = r+offset_r; 

if  (F(p) . BoundaryList (b_cnt , 1)  ==  F(p) . BoundaryList (1 , 1) ) &... 
(F(p) . BoundaryList (b_cnt , 2)  ==  F(p) . BoundaryList ( 1 , 2) ) 
check_cnt  = b_cnt; 

end 

%Get  the  point  before  the  hit 

[back_of f _r ,back_of f _c]  = Neigh8(index-1) ; 

Last_non_boundary  = [r+back_of f _r , c+  back_off_c] ; 

’/.set  up  next  boundary  search  center 
r = F(p) . BoundaryList (b_cnt , 2) ; 
c = F (p) . BoundaryList (b_cnt , 1) ; 
break;  /(exit  the  local  boundary  search  loop 

end 

end 

end 

°/0Delete  duplicate  points 

F(p) . BoundaryList  = unique (F(p) . BoundaryList rows ’ ) ; 

end 

/(To  display  the  merged  regions  first  form  a label  matrix  with  zeros  at 
/(the  feature  boundaries 
NewL_bdr  = zeros (Nr , Nc) ; 

NewL_bdr  = NewL; 
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for  p = 1 :no_f eatures 

ign  = find(ignore  ==  p) ; 
if  “isempty(ign) 
continue ; 

end 

IFpB  = length(F(p) . BoundaryList (:, 1) ) ; 
for  q = l:lFpB 

r = F(p) .BoundaryList (q, 2) ; 
c = F(p) . BoundaryList (q, 1) ; 

NewL_bdr(r ,c)  = 0; 

end 

end 

function  Postprocessor (12 , NewL_bdr) 

yoP0ST_PR0CESSDR  - post-process  feature  boundary  label  matrix 
'/.  Overlay  boundaries  on  cell  image 

% 

% Input : 

'/.  12  - image  pixels 

'/.  NewL_bdr  - feature  label  matrix  with  boundaries  identified  by  zeros 
‘/.Output : 

% Figure  with  black  feature  boundaries  overlaying  the  original  image 

1 

'/.Author : 

'/.  David  E.  Gilsinn 

'/.  Mathematical  and  Computational  Sciences  Division 

'/.  National  Institute  of  Standards  and  Technology 
'/,  100  Bureau  Drive,  Stop  8910 

'/.  Gaithersburg,  MD  20899-8910 

7. 

bdrs  = (NewL_bdr  ==  0) ; 

I2new  = 12; 

I2new(bdrs)  = 1; 
figure , imshow(I2new) ; 
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